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ABSTRACT 

Numerical simulations play an important role in the study of the structure formation of the universe. 
However, the mass resolution in current simulations is still poor. Due to technical difhculties, it is 
necessary to use both a greatly reduced number density of particles and a greatly raised unit particle 
mass. Consequently, the particle masses used in cosmological simulations are about 10^*^ times larger 
than the Gev candidates in particle physics. This is a huge physical bias that cannot be neglected 
when interpreting the results of the simulations. Here we discuss how such a bias affects cold dark 
matter (CDM) cosmological simulations. We find that the small-scale properties of the CDM particle 
system are changed in two aspects. (1) An upper limit is imposed on the spatial resolution of the 
simulation results. (2) Most importantly, an unexpected short mean free path is produced, and the 
corresponding two-body scattering cross section is close to the value expected in the self-interaction 
dark matter model. Since the mean free path of real CDM particle systems is much longer than that 
in the simulations, our results imply that (1) there is probably no "cusp problem" in real CDM halo, 
and (2) a much longer time is needed to form new virialized halos in real CDM particle systems than 
in the simulations. The last result can help us understand the 'substructure problem'. Our discussion 
can also explain why the massive halos in the simulations may have smaller concentration coefficients. 
Subject headings: methods: N-body simulations — cosmology: dark matter — galaxies: halos — 
gravitation 



1. INTRODUCTION. 

Thanks to the rapid development of numerical com- 
putation technology, cosmological iV-body simulations 
have become an increasingly important tool in our un- 
derstanding of the gravitational clustering process in 
the universe. Suto (2003) noted that the number of 
particles used in simulations has risen exponentially 
with time, from about 400 in 1975 to about 10^ to- 
day. This means that current iV-body simulations are 
available to give reliable results on some complex cos- 
mological issues such as large-scale structure formation 
(Einasto,2001) and the properties of dark matter ha- 
los (Coorayfe Sheth,2002;Navarro,Frenk& White,1996). 
Sometimes such numerical methods are the only way to 
get a qualitative picture of the system for which there is 
no mature theory. This is especially the case when strong 
nonlinear evolution on small scales is involved, such as in 
the central region of dark matter halos (Primack,2003) 
and in the processes of galaxy formation and merging 
(Kauffmann,White& Guiderdoni,1993), etc. 

High-resolution simulations based on the ACDM cos- 
mology model have shown nice agreement with the 
large-scale structure observations in several cases (Freed- 
man,W.L.& Turner,M.S. 2003). The predicted power 
spectrum in the simulations agrees well with recent red- 
shift surveys like the Sloan Digital Sky Survey (Tegmark 
et al.,2003) and the Two-Degree Field Galaxy Redshift 
Survey (Percival,W.J. 2001). But the simulations based 
on cold dark matter(CDM) models do not work so well 
on small scales. One aspect is that the central density 
profiles of halos predicted in these simulations are too 
cuspy (Reed D,et al. 2003) to be consistent with the ob- 
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served results of dwarf galaxies and low surface bright- 
ness galaxies (Primack,2003) or to be consistent with the 
X-ray and gravitational lensing study of the distribu- 
tion of dark matter (DM) in galaxy clusters (Haruyoshi 
& Kiyoshi 2004; Tyson et al. 1998). Another as- 
pect is that 10 times more substructures are predicted 
than observed (Willman et al, 2002); this is called the 
" substructure problem" or " dwarf galaxy problem" (Pri- 
mack,2003;Klypin,Kravtsov& Valenzuela,1999). Some 
believe that these problems come from deficiencies in 
the CDM model and therefore introduce some alternative 
models, such as the self-interaction dark matter(SIDM) 
model, warm dark matter (WDM) model, etc. (Martinfe 
Jesper,2003;Celine,et al.,2002;Spergel&Steinhardt,2000) 
Here we find that studying a common difficulty in cur- 
rent TV-body simulations will help us understand the 
above problems. Due to the technical difficulties, the 
total number of particles N we can use even in recent 
simulations is still somewhat limited. The most powerful 
'high-resolution' simulation today with the largest total 
number contains only about 512'^10^ particles. But at 
the same time, we have to keep the mean density of the 
universe at a preset level. So the mass of each particle 
in these simulations has to be as high as about 10^°Mq 
(for a large-scale simulation), comparable to the mass 
of a small galaxy. Now, when we study the large-scale 
structure formation process, it is an acceptable assump- 
tion that each particle in the simulation represents one 
galaxy. However, when we study the small-scale proper- 
ties of dark matter particle systems, such as the CDM 
halo density profile, the merging process of halos and the 
formation of galaxies, the limitations arising from an in- 
adequate N will emerge. For example, the number of 
particles in one CDM halo (^ 10^'^ Af©) in current simu- 
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lations ranges from 10^ to 10^, and so the mass of each 
particle is greater than 10^' Mq, therefore, what we can 
learn from such simulations will concern halos made up 
of such "huge mass particles". 

However, in the view of particle physicists, it is hard 
to accept that the mass of a dark matter particle is some 
10^''' times larger than that of usual GeV candidates. 
And it is also hard for astrophysicists to accept, because 
a point particle with such a huge mass is somehow a mas- 
sive black hole, which will affect the galaxy seriously. If 
we do not believe the mass of the CDM particle should 
be the size of a small galaxy, then the particles in the 
simulations are something different from real CDM par- 
ticles. In other words, we have unavoidably introduced 
a physical bias in current simulations by having to si- 
multaneously reduce the particle number density n and 
raise the mass of each particle m by a huge factor fAE ■ 
The bias is like trying to study the clustering properties 
of a huge society of ants by simulating a colony of ele- 
phants: we call it the "ant-elephant (AE) bias". Because 
CDM particle systems with AE bias can have different 
dynamical and statistical properties, the clustering and 
evolution scenario of CDM halos that we get from the 
simulations may be very different from the truth. One 
might ask, since such a physical bias is present in the 
current cosmological simulations, are the simulation re- 
sults reliable? To what extent will these results reflect 
the actual situation? Is there any connection between 
the bias and the "cusp problem" and the "substructure 
problem" mentioned above? Learning how the AE bias 
affects the simulation results and identifying the circum- 
stances when the effect is serious can help us reconstruct 
a correct scenario of dark matter clustering and evolution 
from the simulation results. 

When the factor Jae is less than about 10^, it is pos- 
sible to test the effects of the AE bias in the numerical 
simulation directly. The number of particles used has 
already increased about 10^ times in the past 30 years, 
which means the factor /ae has already decreased by a 
large amount, yet no serious changes in large-scale prop- 
erties has been reported. The factor Jae can be as big as 
10^*^ in the simulations, and therefore it is indeed neces- 
sary to study how the simulation results will be affected 
by such a physical bias. 

2. PHYSICAL DIFFERENCE CAUSED BY THE AE BIAS 

In this Letter, we study the physical difference the AE 
bias makes on the cold dark matter systems, and find 
that the simulation results are affected in two aspects. 

2.1. Constraint on the Spatial Resolution of Simulation 

Results 

It is easy to sec that the AE bias can impose an upper 
limit Lrnin ou tlic spatial resolution of the simulation re- 
sults. The reason is that the AE bias has greatly reduced 
the particle number density. To get reliable statistical re- 
sults from simulation data and in order to include enough 
simulations particles for consideration the size of the sta- 
tistical sample boxes cannot be; too small. For example, 
when discussing the galaxy formation process at the cen- 
ter region of a CDM halo, one needs to calculate the local 
density distribution p of the CDM particles. Yet, when 
estimating the local density at a certain position, the size 
of the sample box L can not be too small. An very small 



sample box will include very limited simulation particles, 
and the statistical result will be unreliable. 

Consider an extreme case: if we use an ultra-small sam- 
ple box with volume Vi ^ nisim/p, then we can only find 
about one simulation particle in it. Here the one simula- 
tion particle with huge mass nisim actually represents all 
the real CDM particles within the sample box. This is 
not a good approximation but it is unavoidable because 
of technical difficulties, i.e. we are unable to deal with 
so many particles. Consequently, the simulation results 
can-not provide reliable information when the physical 

processes occur on scales smaller than Lmini ~ V^^^ . 
We have to be aware of this upper limit L„iini when us- 
ing the simulation results to discuss the galaxy formation 
process and the merging process of substructures. 

Lmini can be estimated in this way: the simulation 
particle number within a sample box with size L is about 
N ~ pL^/rUsim, where nisim is the mass of each simula- 
tion particle and p is the local density; set A'' = 1 and we 
get 

\ p J 

Here we represent msim and p with the solar mass( 
Mq) and the critical density of the universe pcr = 
■iH^/8nG = 1.88 x lO^^e . . Equation(l) gives 
a physical limitation on the spacial resolution of the sim- 
ulation results. 

However, it is still not enough, for it is hard to do 
statistical analysis using a sample box that can include 
only about one particle. In general, when we discuss 
the statistical results of a given simulation sample, the 
statistical fluctuation must be taken into account. For 
example, when we calculate the local density p at the 
halo center from the simulation result, the volume of the 
sample box must be big enough to be able to contain 
adequate simulation particles. The least particle num- 
ber that one sample box must contain depends on the 
requirement of a reliable statistical result. A sample box 
with a size round Lmini is still not big enough and it can 
hardly give a reliable statistical result of p. 

Let us consider the statistical result of p being pro- 
portional to the particle number found in the sample 
box with given volume V: then the statistical fluctua- 
tion of p is about Sp/p cx SN/N. For a Gauss random 
field, SN/N ~ 1/y/N. Suggest a reliable statistical re- 
sult requires Sp/p cx SN/N ~ 1/Vn < a%, where N is 
the particle number found in the sample box with given 
volume V: then we get N > 10^/a^. This result means 
that the sample box must include at least N > 
particles in order to give a reliable statistical result of p. 
Thus, we get another limitation on the size of the sample 
box L„,„2 = 1/1/3 > (A^Fi)! ~ L„„„i(10Va')^.^ 

Lmin2 givcs a Statistical limitation on the spatial reso- 
lution of the simulation results. Here the factor a varies 
between and 100 based on different problems: Lmin2 
is commonly larger than L^nini- Only when a — > 100 
is imm2 ~^ ^mini- On scalcs Smaller than Lmin2, not 
only can the statistical error be serious but the phys- 
ical process can be affected as well. Let us take the 
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study of the galaxy formation process at the halo cen- 
ter region as an example again; one can find few massive 
particles in simulation using a sample box smaller than 
Lmin2, but One can find numerous real CDM particles 
that can be smoothly distributed in this field. There- 
fore the simulation of the CDM particle system will of- 
fer much more asymmetry and instability to the gravi- 
tational background for the study the galaxy formation 
process on such a scale. 

Lmini and Lmin2 Can help us estimate a lower bound 
to the scale of applicability of the simulation results. For 
example, let us suggest that we have a ACDM cosmology 
simulation with rrisim = IO^Mq. For the center region 
of one CDM halo, p/pcr ~ 10 , set a = 1 and we get 
Lmini ~ Ikpc/h, Lmin2 > lOkpc/h. When using these 
simulation results to discuss the central properties of the 
halo, such as the density profile or the galaxy formation 
process, the AE bias must be taken into account. This 
is not an exciting conclusion, and we have to accept the 
fact that by raising the mass resolution 10'^ times, the 
spatial resolution will only increase by a factor of 10. 

This AE bias effect comes from the statistical error 
SN/N, which only occurs on small scales and can be 
neglected on large scales where plenty of simulation par- 
ticles are considered to reach a reliable statistical result. 

2.2. Variation of the Mean Free Path 

The second effect of the AE bias is a variation in the 
mean free path. Since gravitation is a long-range inter- 
action, there are three different factors that can change 
the CDM particle velocity. One is the gravitation of the 
existing halo, i.e., after the halo has formed. The mean 
free path of the particles L will then be about the size of 
the halo, R. Another factor is the scattering of two CDM 
particles. Suppose that particle 1 passes by a stable par- 
ticle 2 with velocity v. If particle 1 gains a scattering 
angle larger than 45° (in this case \Sv\ ~ and we call 
it a 'scattering'), the impact parameter b should satisfy 
b < 2Gm/v^, where m is the particle mass. Let n be the 
particle number density and set the equal to the veloc- 
ity dispersion cr^; then the mean free path of scattering 
Ls should satisfy n(Ls7r6^) ^ 1, and we get 

Ls = at/{i7TG^mp) (2) 

If we define T^cr^ — Lg, we can also expect the two-body 
relaxation time T„ to be proportional to a^/G^mp oc 
1/m (X N (for a given and density p). The third 
resulting proportionality comes from summing the ve- 
locity deflections from distant particles. Binney & 
Tremaine(1987) estimated the relaxation time that was 
due to this cause to be proportional to N/lnN. But 
analytical and simulation results (Huang et al.l993; Die- 
mand et al.2003) show that the relaxation time is pro- 
portional to N, not N/lnN, if the softness factor is taken 
into account in the simulations. Here we use eq l(2Jl to 
discuss the effect of the AE bias. 

As we can see from the formula above, for certain val- 
ues of p and CTu , the AE bias will raise the value of m and 
make Lg shorter. For a virialized halo, we can estimate 
the velocity dispersion from galaxy dynamics (Binney 
&Tremaine,1987); cr^ ~ GM/R, with M the total mass, 
and R the radius of the halo. Substitute a„ into eq.©, 
and we have 
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where p — M/jTtR^. If the AE bias causes m — > mfAE 
and N N/f^E, then we get 

Ls Ls/ Jae and — > Ty/ Jae (4) 
For example, in the central region of a typical halo in 
the simulation, let R = IMpc/h, M/m — Nhaio = 10^, 
and p/p ~ 10^, and we then get Lg = 3.3Mpc/h, close 
to the dimension of the halo. If we apply the expres- 
sion of the scattering cross section commonly used in the 
SIDM model, asi = l/{Lp), and if we let the central 
area density of the halo be p « lOOGev/cm^, we can 
estimate aAE ~ 9 x 1Q~^^ cm^ /GeV. In the above, we 
only used typical data of the central area of halos. In- 
deed, the mean scattering cross section of the whole halo 
should be smaller. However, considering that the mass 
of the DM particle candidates is on the order of GeV, 
/ae ^ 10™, from eq.Q we can see that the actual scat- 
tering cross section between two DM particles should be 
far less than the cross section in the simulation. In this 
case, the change of the particle velocity depends mostly 
on the deflection by the existing halo. When discussing 
the SIDM model, Spergel &Steinhardt(2000) and Dave 
et al.(2001) suggested that if the asi of the DM particles 
is around 10"^'^ to 10~^'^cm^GeV~^ , or their mean free 
path L ^ Ikpc — IMpc is at the solar radius level, then 
the halo center cusp problem and the excessive substruc- 
ture problem can be solved. Now we can see that the AE 
bias has provided an unexpected short mean free path 
Lae and a ctae ^ crs/ smaller than the SIDM model 
expectation. 

To give a rough estimate here, we set the simulation 
timescale to be taii ~ lOOGyr, and the particle velocity 
V ~ lOOOfcms"^, and Lg ^ R ^ IMpc. Then one par- 
ticle in simulation will scatter for about taii/tgcatter — 
tall/ {Lg/v) ~ 10 times. Yet the real CDM particles will 
not experience such scattering for they have a far more 
smaller scattering cross section. This is the difference 
caused by the AE bias. 

Unfortunately, at present, there is not a quantitative 
theory on how the mean free path Lg will affect the his- 
tory of structure formation, but the qualitative picture is 
clear. When Lg R, a short Lg means a large scatter- 
ing cross section and a high thermal conductivity, which 
might alleviate the cusp problem and make the halos 
more spherical. But some authors argue that a short Lg 
can also drain energy from the halo centers, which could 
lead to core collapse (Burkert,2001). When Lg ^ R, 
a large Lg means a small scattering cross section and a 
large average relaxation time T„. Eq.Q shows us that 
with improvement in the mass resolution (decrease in 
/ae), Tv will be increased. The simulation by Diemand, 
et al.(2004) shows that Ty oc N and that it can easily rise 
to more than 1000 Gyr by increasing the total particle 
number TV in the simulation. In this case, the tidal force 
is unable to damage small halos in one Hubble time, as 
was shown in the simulation by Kazantzidis,et al.(2003). 
Here the substructures are mainly what remain of halo 
merging. A very small scattering cross section also means 
that the effects of the two-body scattering of the CDM 
particles can be neglected and that a very poor ther- 
mal conductivity is expected. If is far much longer 
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than one Hubble time, then it is hard to make the CDM 
halo virialize with the poor thermal conductivity. And a 
new density distribution implies a different particle en- 
ergy distribution, so there will not be enough time to 
form new substructures on scales less than R, such as 
the p ~ cusp profile in the isothermal model. 

The current simulations have brought in an unexpect- 
edly short Ls ~ i?, not in both cases above, and it is 
different from the real CDM particle systems (Lg ^ R). 
The corresponding scattering cross section is smaller 
than what is needed in the SIDM model, and it can not 
be neglected. The simulation results will be sensitive to 
Ls and to the relaxation effects of subhalos, such as the 
number and mass density of subhalos (Ma et al.2004). 
One can expect that, with improvements to the numer- 
ical technique, the AE bias will be gradually alleviated, 
and the simulation results will gradually move to the 
Ls R regime. By increasing Ls, the halo center cusp 
problem can be alleviated. But the relaxation time will 
soon rise up to more than one Hubble time. As a result, 
halos will tend to be more triaxial, and more small ha- 
los will survive the process of merging. In other words, 
the AE bias decreases the number of substructures. The 
substructure problem in future simulations can become 
more serious. Actually, this is in agreement to a certain 
extent with the discussion of gravitational lensing (Mao 
et al.2004) in which more substructures are expected in 
real DM halos than in the current simulations. But if 
Ls increases more, the thermal conductivity in the halos 
will become poor, which means it will be harder to form 
enough virialized halos in one Hubble time. 

By the way, eq.Q can help us understand the con- 
centration properties of DM halos in simulations. In a 
given simulation in which the particles have a fixed mass 
m, bigger halos will have larger total masses M. Then, 
based on eq.JHJ, one can predict that the big halos will 
have larger Lg/R than small halos. Large Ls/R mean 
that the particles can move more freely in the halo, which 
will become more dispersed, so equation(3) can give a ra- 
tional explanation of why the more massive halos in the 
simulations have a smaller concentration coefficient C . 
Eq. is available for all virialized particle systems with 
a pure Newton gravitation interaction. If we suppose 
the particles to be the stars of elliptical galaxies, then 
the discussion above will predict that the more massive 
elliptical galaxies will be more dispersed. In other words, 
they will have a fainter mean surface brightness, and this 
is consistent with the observations (Binney &Merrifield, 
1998,p.205). 

3. CONCLUSION AND DISCUSSION 

In summary, we have pointed out that a physical bias 
exists in current cosmological simulations. We find that 



Binney, J. &Merrifield, M. 1998, GaZaciic Astronomy, ( Princeton 

University Press, P205 ) 
Binney, J. &Tremaine,S. 1987, Ga/aciic Dynamics, { Princeton 

University Press ) 
Burkert,A.2001,to appear in DARK 2000 

(Springer) astro-ph/0012178 
Celine,B. et al.,2002,PHD, 66,3505 
Cooray,A.&Sheth,R. 2002,PhR,372,l 
Dave,R. et al.2001,ApJ, 547,574 
Diemand.J. etal. 2004,MNRAS,348,977 



the AE bias can change the small-scale properties of the 
CDM particle system in two respects. (1) An upper limit 
is imposed on the spatial resolution of the simulation re- 
sults. (2) An unexpected short mean free path of the sys- 
tem is generated. The corresponding two-body scatter- 
ing cross section is near the value expected in the SIDM 
model. 

What can we learn from the discussion above? The 
upper limit imin2 of the first effect can help us estimate 
the reliable range of the simulation results. The second 
effect can alter the simulation results on small scales and 
cause the halo center to be more cuspy and more spher- 
ical. These results can help us understand the "cusp 
problem" and the " substructure problem" in the current 
simulations, and they can help us reconstruct a correct 
scenario for the clustering and evolution of the real dark 
matter in the universe. Because the AE bias factor Jae 
is currently about 10™, the mean free path Ls and the 
relaxation time of real CDM particles systems should 
be much larger than we see in the simulations. When Ls 
is far more larger than the halo diameter, the halo will 
be dispersed and will be hard put to form a cuspy den- 
sity profile. When T„ is much longer than one Hubble 
time, CDM halos will not have enough time to be viri- 
alized of the poor thermal conductivity (this is contrary 
to the prevailing view of the popular halo model), and 
nearly all the centers of the existing halos can survive the 
process of merging until the present time. Our results 
imply (1) that real CDM halos may not have the cusp 
problem(the cusps of the simulation halos are mainly the 
effect of the AE bias)and (2) that real CDM halos may 
not have enough time to be virialized. This result can 
help us understand the substructure problem. The dis- 
cussion above also tells us that the theoretical analysis 
of the halo property, which is based on plenty of thermal 
conductivity, is available for simulation halos but not for 
the real CDM halos. 

On the other hand, since it is hard to form new viri- 
alized halos in one Hubble time and since most existing 
small halos (coming from primary disturbances) can sur- 
vive the process of merging, the number of these small 
halos may retain some information about the primordial 
disturbance. If we believe the number of dwarf galax- 
ies is proportional to the number of small halos, then 
the percentage of dwarf galaxies in galaxy clusters can 
tell us something about the small-scale properties of the 
primordial disturbance. 
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